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The linear stability of rapid granular flow on a slope under gravity against the longitu- 
dinal perturbation is analyzed using hydrodynamic equations. It is demonstrated that 
the steady flow uniform along the flow direction becomes unstable against the long- 
wavelength perturbations longitudinal to the flow direction for certain parameter ranges 
to form the density wave, in contrast with the finite wavelength instability against the 
transverse perturbation (Forterre & Pouliquen 2002) . It is shown that the instability can 
be understood as the the long-wave instability of the kinematic waves in a quasi-one 
dimensional system. The results are compared with our previous molecular dynamics 
simulations (Mitarai & Nakanishi 2001), where the spontaneous density wave formation 
has been found. 



1. Introduction 

Granular flow exhibits a variety of dynamical phenomena, which have been attract- 
ing research interests for years (for reviews, see e.g. Savage 1984 and Jaeger, Nagel & 
Behringer 1996). Its complex behaviors can be seen even in a simple situation like the 
gravitational flow on a slope. When the inclination angle is large and the slope is rough, 
a rapid and relatively low-density flow is realized, and the interaction between grains is 
dominated by inelastic collisions. On the other hand, when the inclination angle is small, 
the flow becomes dense and slow, and the frictional interaction plays an important role 
(Savage 1984; Mitarai & Nakanishi 2003). The comprehensive rheology of the granular 
flow has not been established yet, except for the rapid collisional flow regime, where the 
hydrodynamic models have been developed with the constitutive relations based on the 
kinetic theory of inelastic hard spheres (Jenkins & Savage 1983; Campbell 1990; Lun 
et 0,1.1984; Goldhirsh 2003); it has been demonstrated that a certain quantitative agree- 
ment can be achieved for the steady flow by introducing the spinning motion of each 
grains (Mitarai, Hayakawa & Nakanishi 2002). The steady granular flow, however, turns 
out to be unstable in various ways, and shows rich phenomena. 

In the experiment on a shallow granular flow on a wide slope, [Forterre fc Pouliquen(2001) 
have observed that there appears a regular pattern of longitudinal streaks along the 
flowing direction. This phenomenon has been analyzed by means of the hydrodynamic 
equations for rapid granular flow (Forterre & Pouliquen 2002). They have calculated the 
steady solutions and examined its linear stability numerically. They have found that, 
at a certain parameter region, the steady flow shows the "inverted density profiles" , in 
which the maximum density appears not at the bottom but at a finite distance from the 
bottom because of the agitation by the collisions with the rough solid bottom. It has 
been shown that the solutions with the "inverted density profiles" are unstable against 
the perturbations transverse to the flowing direction, and the instability results in the 
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vortex patterns analogous to the rolls in the Raylcigh-Bcrnard convection; the streaks 
found in the experiment were interpreted as the result of the rolls of vortices. 

Another instability that has been observed is the density wave formation along the 
flowing direction; experimentally, Louge & Keast(2001) have observed the jamming pat- 
ters traveling upstream in the dense chute flow, and Prasad, Pal & Romkens (2000) have 
found the waves develops in the shallow flow as they travel downstream. The present 
authors have performed the molecular dynamics simulations of two-dimensional granular 
flow on a slope and found that the steady flows are unstable against the density wave 
formation when the length of the slope is long enough and/or the particle density is low 
enough (Mitarai & Nakanishi 2001) . 

The purpose of this paper is to perform the linear stability analysis on the hydro- 
dynamic equations to investigate the nature of the density wave formation instability 
found in the experiments and the numerical simulations. The basic method is the same 
with that used in Forterre fc Pouliquen(200"2){ but we examine the stability against 
the perturbations longitudinal to the flowing direction in contrast with the case of 



Forterre & Pouliquen(2002) where the transverse stability has been studied. 



This paper is organized as follows. In the hydrodynamic model for rapid granular 
flow is introduced. The steady solutions uniform along the flow direction are numerically 
obtained in fJ3] and the results of the linear stability analysis is presented in 2J In 
discussions and the comparison with the molecular dynamics simulations are given. The 
results are summarized in Sj^l 



u ■ VJ p = -pS7 ■ u, (2.1) 
u ■ V ) u = pg - V • E, (2.2) 



2. Hydrodynamic Equations for Granular Flows 

2.1. Hydrodynamic equations and constitutive relations 

The hydrodynamic fields for granular flows in three dimensions are the mass density p, 
the mean velocity u, and the granular temperature T, where T =< Su 2 > /3. Here, 
5u = u— < u > and < . . . > represents the average over the microscopic scale. Under 
gravity, they follow 

d_ 
dt 
d_ 

dt ' 

^p(J i +u-V S )T=-V-q-X:Vu-T ) (2.3) 

with the acceleration of gravity g, the stress tensor E, the heat flux q, and the energy 
loss r due to the inelastic nature of interactions between grains. 

We employ the constitutive relations derived by Lun et aL(1984) for three dimensional 
system based on the kinetic theory of the inelastic particles: f 

E = (p - CV • u)I - 2pS, (2.4) 
q = -kVT, (2.5) 

f The original form of q derived by Lun et aZ.(1984) is q — —kVT — k^Sv. The coefficient 
Kh is proportional to (1 — e p ), thus disappears in the elastic limit. We checked that the influence 
of the term KhSlv on t he steady solutions is small in the considered parameter region, therefore 
neglected this term as |Forterre fc Pouliquen(2002)| 
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Table 1. Dimensionless functions used in the constitutive relations and the boundary 

conditions with rj = (1 + e p )/2. 



where 

S= -[Vw + (Vu)*] - -(V-m)I, (2.6) 

and 

p(i/, T) = p p A(^)T, T) = Pp af 2 ( V )T^ 2 , T) = p p af 3 (u)T^ 2 , 

K{v,T)=p p o-U{v)T 1 '*, T{ V ,T) = ^{\-el)h( V )T^\ (2.7) 

a 

with the material density of particle p p , the packing fraction v = p/p p , the particle 
diameter a, and the restitution coefficient between particles e p . Here, / represents the 
unit matrix. The dimensionless functions fi(y) (i = 1, ... ,5) are given in Tabled F° r 
the radial distribution function g${v) in these functions, we adopted the form suggested 



by Lun fc Savage(1986) 



goH = (i - *K0»- ' (2 - 8) 



with the maximum solid fraction v m , for which we use 0.60 as Forterre & Pouliquen(2002) 
In the following, all the variables are non-dimensionalized by the length unit a, the 
mass unit p p o- 3 , and the time unit yfajg. The density field is expressed by the packing 
fraction v instead of the mass density p. The restitution coefficient between particles is 
set to be e p = 0.7, that is the value used in our previous simulations (Mitarai & Nakanishi 
2001). 

2.2. Boundary conditions 

The granular flow has a non-zero slip velocity at the solid boundary, where we should 
impose the momentum and the kinetic energy balances. Sophisticated boundary con- 
ditions have been proposed based on microscopic calculations of the kinetic theory for 
specific geometries (Jenkins & Richman 1986; Richman 1988; Jenkins 1992). We employ, 
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however, a simpler form of the boundary condition obtained from a heuristic approach 
(Johnson & Jackson 1987; Johnson, Nott & Jackson 1987); 

-n-£.t = »j*(i/,T)|u.|, (2.9) 
n-q = -u a -E-n-T*(v,T), (2.10) 

where the unit vector n is normal to the floor, u s is the slip velocity at the floor, and 
t = u s /\u s \ is the unit vector along the slip velocity. 

The first equation l|2.9|l expresses that the stress at the boundary balances with the 
momentum transfer due to the collisions between the slope and the flowing grains. The 
momentum transfer, or RHS of l|2.9|) . is assumed to be given by 

» 7 > ) T)|u.| = ^|ti a |n(i/,T), (2.11) 

with the collision rate Q(v, T) per unit time per unit area. Here, the factor 7r/6 comes 
from the non-dimensionalization of the particle mass m = <t 3 7t/6. The parameter <f> 
characterizes the roughness of the boundary, and the expression means that the fraction 
4> of the particle momentum is transfered to the boundary in each collision, therefore, 
the larger value of <f> represents the rougher boundary. For the rough boundary in the 
two-dimensional system with closed packed disks, Jenkins & Richman(1986) estimated 



as sa 0.1, but Forterre & Pouliquen(2002) adopted the smaller value <j> = 0.05 for most 



of the cases because they expected that a boundary with closed packed spheres in three- 
dimensional system is smoother on average. In this paper, we mainly use (f> = 0.05, but 
the case of = 0.10 is also examined in order to see general tendencies. 

The second equation i|2.10[l represents the energy balance, and means that the heat 
flux at the boundary comes from two effects, namely, the frictional heating due to the 
non-zero slip velocity and the energy loss due to the inelastic collision with the floor. The 
energy loss term T* in (|2.10|l is given by 

r* = $-|-|r-^,T) (2.12) 

with the collision rate fi(z/, T). The parameter <£> represents the rate of energy loss per 
collisionf , and we use $ = 0.39 in this paper. 

For the collision rate Q(v, T), we use the form adopted by Forterre & Pouliquen(2002) 



i.e., fl(v, T) = yWTvgoiy) / v m : then the expression of rj*(u,T) and T*(u,T) in non- 
dimcnsionalized form are 

rf{v,T) = 4>f 7 {v)h{v)T 1 /\ (2.13) 

r>,T) = i$/ 6 H/ 4 HT 3 / 2 , (2.14) 

where the dimensionless functions fe{v) and fi{v) are given in Tabled 

At infinity, we impose the condition that the stress and the heat flux vanish, namely, 

S^0 and q — > as y — > oo, (2-15) 

where the y axis is taken perpendicular to the floor (figure^. 

| |Johnson fc Jackson(1987)1 explicitly relate $ to the restitution coefficient between the floor 
and the particles e w in the lorm $ = (1 — ejj,), but we adopt 12.12H as a more general expression. 
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Figure 1 . The coordinate system. 



3. Steady Flows 

3.1. Equations and Numerical method 

First, we consider the steady solution which is uniform along the slope for (|2.1() - 12. .'{1 
with the boundary conditions (|2.9|l . (|2.10(l . and Q2.15D in the form 



v{x,y,z,i) = v (y), 
u(x,y,z,t) = {u (y), 0,0), 
T(x,y,z,t) = T (y), 



(3.1) 
(3.2) 
(3.3) 



where we take the x axis along the slope, the y axis perpendicular to the floor, and the 
z axis perpendicular to the x — y plane (figure ^) . 
Then, the equations 12.1fl - l|2.3[) are written as 



= isq sin ( 



dy 

a2 ^yy 



= 



-vq cos 6 

dy 

du dgg 
xy dy dy 



(3.4) 
(3.5) 
(3.6) 



where the superscript denotes that the functions are for the steady solution, namely, 
q a = — k{vq, To)VT = (0, g°, 0), etc. By integrating (|3.4|) and (|3.5|1 over y with the stress 
free condition at infinity, we obtain the condition 



^xy — ~ tan# E° y . 
From 13.4|) - (|3.7|) and the constitutive relations, we have 

f°T ' + is cos6 



v' (y) 



u' (y) = ^^tauh». 



fS 



1 



TO = to 



I 1 e p)J5 1 J2 U J^v v 1 Ji 



2T 



(3.7) 

(3.8) 

(3.9) 
(3.10) 



where /? = fi(uo), f° v = ^fi{ v )\v=v i and the prime indicates the derivative by its 
argument. 

The boundary conditions (|2.9|) and (|2.10|) at the floor (y = 0) for the steady solution 
can be written as 



T 



I/O 



(3.11) 




Figure 2. (a)The contour lines of 6 in the (h, v) plane for <f> — 0.05. The region of the 
non-monotonic density profiles is shown by a grey region. (b)The flux Qo vs the one-dimensional 
density po for (j> — 0.05. The solid line and the dashed line are for 6 — 14° and 8 — 16°, respec- 
tively. 



K = -/e° ~ \® T <) ■ (3-12) 

The boundary condition 12.15fl that the stress and the energy flux should vanish at 
infinity is satisfied when (Ahn, Brennen & Sabersky 1992) 

T^(y)~>Q when y -> oo. (3.13) 

In order to obtain the steady solutions, we integrate H3.8(l . (|3.9[) . and H3.1()(l numerically 
using the fourth-order Runge-Kutta method with the boundary conditions . 1 If) and 
<|3.12[1 . We employ the shooting method to find the solution which satisfies the condition 
(|3.13|l (Forterre & Pouliquen 2002); for a given inclination angle 9 and a given density 
at the floor z^o(0), we search for a solution by adjusting the value of the velocity at the 
floor Uo(0). In the actual calculations, we integrate the equations numerically from y = 
to a certain height y max , and search for the solution which gives |TQ(y max )| < 10~ 7 . The 
value of j/max is chosen to be large enough in comparison with the relaxation length, that 
depends on the parameters and can be determined only after the solution is obtained. 

We use 9 and i'o(O) to specify the solution in the rest of the paper. 

3.2. Numerical solutions 

For a given roughness <f> of the slope, steady solutions are found for a certain range of 
the inclination angle of the slope 9 (Forterre & Pouliquen 2002). We present the steady 
solutions for the two cases, (i) 0=0.05 and (ii) 4> — 0.10; the most of the results are for 
the case of <f> — 0.05, and the case of <fi = 0.10 will be given to examine general trend. 

3.2.1. The case of <f> = 0.05 

For an appropriate 9, there exist steady solutions for a given density at the floor 
^o(0). We numerically find the steady solution for the range 8° < 9 < 20° for the 
moderate density; towards the lower limit of 9, the length scale of the density decay in 
the y-direction tends to zero and becomes smaller than the particle diameter, which is 
physically unacceptable, while the decay length of density diverges towards the upper 



limit of 9. This is consistent with the analysis by Anderson & Jackson(1992) that the 
steady solution in the high-density limit is allowed for a finite range of 9. 

For a given steady solution, we define the "thickness" h and the "mean density" v: the 
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Figure 3. The steady solutions for <f) — 0.05 and 6 — 16°; the density profiles (a), the velocity 
profiles (b), and the temperature profiles (c). Different line types correspond to the solutions 
with different densities; ^o(0) = 0.04 (solid lines), 0.05 (long-dashed lines), 0.08 (short-dashed 
lines), O.fO (dotted lines), and 0.f5 (dash-dotted lines). 



thickness h is the value of y where the density is 1% of the maximum density, and the 
mean density v is given by 

1 /*°° 

v (y)dy. (3.14) 



It has been found numerically by Forterre & Pouliquen(2002) that there exists a one to 



one correspondence between (9,uq(0)) and {h,v). 

Figure[5Ja) shows the contour lines of 9 in the (h, v) plane, where 9 increases from left 
to right. In this plot, it is clear that h goes to zero as 9 decreases and h diverges as 9 
increases. There is a separatrix near the bottom at h « 7, where the value of 9 is around 
15°; the contour line for 9 ^ 14° leads to h = as v decreases, while h diverges along 
the contour lines for 9 ^ 15° as v becomes small. 

The difference between 9 ^ 15° and 9 ^ 14° can be seen in the relation between the 
flux Qq and the one-dimensional density po defined as 

/>oo />oo 

Qo = / v (y)u {y)&y and po = v Q {y)Ay, (3.15) 
Jo Jo 

respectively, for a fixed inclination angle. It is found that Qo is an increasing function of 
po for 9 ^ 14°, while Qq has a minimum at a finite po for 9 ^ 15°; the plot of Qo vs po 
is shown in figure Hfb) for 9 = 14° and 9 — 16°. 

The typical profiles of the density, the velocity, and the temperature for 9 = 16° arc 
shown in figureOlfor the density at the floor ^o(0) = 0.04 ~ 0.15. We see in figure |^a) 
that the density decays monotonically when the density at the floor ^o(0) is small enough 
(z^o(0) ^ 0.10), while for higher density (^o(0) = 0.15) the maximum density appears at a 
finite height. The region where the maximum density appears at a finite height is shown 
in figure [2t a) as a grey region. We focus on the lower density region because the density 
decays monotonically in our previous simulation (Mitarai & Nakanishi 2001; Mitarai 
et a/.2002). 

For ^o(0) < 0.10, the higher density flow shows lower flowing speed in the case of 
9 = 16°, which results in the decrease of the flux Qo as Po increase for po < 0.2 (figure 
12b)). For higher density (^o(0) = 0.15 in figureOJ), the flowing speed gets faster with p , 
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Figure 4. (a)The contour lines of 6 in the (h, v) plane for <f> = 0.10. The region of non-monotonic 
density profiles is shown by a grey region. (b)Relation between Qo vs po for (f> — 0.10. The solid 
line and the dashed line are for 9 = 20° and 8 = 21°, respectively. 
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Figure 5. The steady solutions for <f> = 0.10 and 9 = 20°. The density profiles (a), the velocity 
profiles (b), and the temperature profiles (c). Different line types correspond to the solutions 
with different densities; fo(0) = 0.10 (solid lines), 0.15 (long-dashed lines) 0.17 (short-dashed 
lines), and 0.20 (dotted lines). 



which causes the increase of the flux Qo for higher po. As a result, Qq has a minimum 
at a finite density. 

In the case of 8 $C 14°, the velocity continuously decrease as the density becomes lower, 
and Qq becomes a increasing function of pq. 

3.2.2. The case of <f> = 0.10 

The slope is rougher than the previous case, and the steady solution exists for 12° < 
9 < 25°. The contour lines for 9 in the (h, v) plane are shown in figure Ufa). As in 
the case of </> = 0.05, h goes to zero as density becomes lower for smaller 9, while h 
diverges for smaller density when 9 is large enough. Figure 0{b) shows p~o dependence 
of Qo-, which is a monotonically increasing function for 9 = 20° and has a minimum for 
9 = 21°. The typical solutions are shown for 9 = 20° in figure [S] For large enough h and 
v, the maximum density appears at a finite distance from the floor. The region of the 
non-monotonic density profile is shown by a grey region in figure Ufa). 
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4. Linear Stability Analysis; Density Wave Formation 

4.1. Normal mode analysis 

We restrict our stability analysis to the perturbation uniform along the z direction, 
because we are interested in the instability along the flow direction. The flow is perturbed 
around the steady solution as 



The governing equations and the boundary conditions are linearized with respect to the 
deviations V\, u\, v\, and T%, but the resulting expressions are rather lengthy and given 
in Appendix A. 

Now we look for the normal modes for the density, the velocity, and the temperature 
perturbations of the form 



The flow is linearly unstable if Re(a) > 0. 

As for the boundary condition at the free surface, the asymptotic behavior of the 
perturbations at large y are used. When y is much larger than the decay length of the 
density and thus uo(y) is very small, the density perturbation should also decay (y cx vq), 
and u, v, and T decay proportional to exp(— ky) (Forterre & Pouliquen 2002), therefore, 
we imposed the boundary condition that 



at y = y max ; fcax is a large enough height where vo(y m &x) < 10 -9 , in addition to the 
condition |To(y max )| < 10~ 7 discussed in 

We solve the eigenvalue problems of the linearized equations (|A 1|I - (|A 4|) for (|4.4|l nu- 
merically using the Chebychev collocation method with the discretization in the y di- 
rection (Gottilieb, Hussaini & Orsag 1984; Canuto et aZ.1988; Boyd 2001; Foterre & 
Pouliquen 2002). It is known that the straightforward discretization of space requires 
two extra boundary conditions (Malik 1990; Foterre & Pouliquen 2002), for which we 
use the momentum balance condition in the y direction at y = and the decay condi- 
tion for the density perturbation, i.e., v'{y) = — (cos 0/To(y))£>(y) at y = y max - In actual 
numerics to solve the generalized eigenvalue problem in the form AV = aBV for the 
complex eigenvalues a and the eigenvectors V, we used LAPACK version 3.0 (Anderson 
et al. 1999). The number of discretization Nd is about 100. 

The numerically obtained eigenmodes contains unphysical modes, called spurious modes, 
due to the discretization (Mayer & Powell 1992; Boyd 2001; Forterre & Pouliquen 2002). 
For the spurious modes, it is known that the Chebychev coefficients of higher wave num- 
ber are large, and the eigenvalues are sensitive to small change of Nd- We determine the 
eigenmodes as physical ones by checking that their Chebychev coefficients for higher wave 
number are small and their eigenvalues varies little upon changing Nd'. We confirmed that, 
for these modes, the highest ten coefficients are less than 10~ 7 when the eigenvectors are 
normalized so that the sum of the absolute values of the real part and the imaginary part 
of the largest component becomes one, and the variation of the eigenvalues by the small 
change of Nd are less than 10 -7 . 



v(x,y,t) = v {y) + ui(x,y,t), 

u(x,y,t) = (u (y), 0,0) + {u 1 (x,y,t),v 1 (x,y,t)),0), 
T(x,y,t) = T (y) + T 1 (x,y,t). 



(4.1) 
(4.2) 
(4.3) 




(4.4) 



u'(y) = -ku(y), v'(y) = -kv(y), f'(y) = -kf(y), 



(4.5) 
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h 

Figure 6. The stability diagram for h vs v for <f> — 0.05, where the unstable (stable) regimes 
are shown by grey (white) regions. We find three unstable regimes A, B, and C. The dashed line 
shows the boundary between the region of the non-monotonic density profiles and that of the 
monotonic density profiles (see also figure The dispersion relations at the points P1~P5 are 
shown in figures |7| ~ |3] The region of h < 1 where the stability boundary is shown by a dotted 
line was not examined in detail because the flow with the decay length less than the particle 
diameter is physically unacceptable. 



(a) (b) (c) 

Re(a) Re(a) Re(a) 




0.01 0.02 0.01 0.02 0.01 0.02 



k k k 

Figure 7. The dispersion relations of the least stable mode for the steady solutions with <f> = 0.05 
and 6 = 16° for i/ (0) = 0.05 (a), 0.08 (b), and 0.10 (c), which correspond to the points PI, P2, 
and P3 in figure |5] respectively. 

4.2. Stability diagram, dispersion relations, and eigenf unctions 

We present the results of the linear stability analysis in the cases of cj> = 0.05 and <f> = 0.10, 
for which the steady solutions are shown in 

4.2.1. The case of <j> = 0.05 

The stability diagram is shown in figure[B]in the parameter space of h vs v. The unstable 
(stable) regimes are shown by grey (white) regions, and the boundary between the region 
of the monotonic density profiles and that of the non-monotonic density profiles is shown 
by a dashed line. Within the investigated region, we find three unstable regimes; the 
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Figure 8. The magnification of the dispersion relations near k = of the unstable stable mode 
for the steady solutions with <j> = 0.05 and 9 = 16° for ^o(0) = 0.05 at PI in the regime B. Re [a] 
grows quadratically in k, while the slope of Im[a] is given by — dQo/dpo in the long-wavelength 
limit, which is shown by a dashed line. 
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Figure 9. The dispersion relations of the unstable mode for <j> — 0.05 with (a) 9 — 14° and 
ID (0) = 0.05 (P4; in regime A) and (b) 9 = 17° and v (0) = 0.10 (P5; in regime C). 



regime A at small h and small v region, the regime B at large h and small v region, and 
the regime C at large h and large v region. When we decrease the density with a constant 
inclination angle (along a contour in figure |2J) , eventually we will encounter either the 
regime A or the regime B, namely, the flow with low enough density is always unstable. 
The two regimes are different in the steady flow behavior as we have seen already: in the 
regime A at the small h side, the flow becomes slower as the density becomes smaller, 
while the flow in the regime B at the large h side flows down faster for the smaller density. 
On the other hand, the denser flow can be unstable in the regime C, which lies within the 
region of the non-monotonic density profiles. In this regime, the denser flow goes faster 
as in the regime A. 

The dispersion relations of the least stable modes a — a(k) are shown in figure [7| for 
= 16° and i/ (0) = 0.05 (a), 0.08 (b), and 0.10 (c), which correspond to the points PI, 
P2, and P3, respectively, in the stability diagram of figure EJ PI lies within the unstable 
regime B. In all cases, it is found that the least stable mode satisfies a(0) = 0. The 
growth rate Re(a) is positive for fo(0) = 0.05 for < k ^ k c with k c w 0.007 (figure 
da)). The magnification around k = in figure El shows that Re(a) grows quadratically 
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Figure 10. The eigenfunctions of the least stable mode for <j> = 0.05, 6 = 16°, and uq(0) = 0.10 
(P3 in figure ^J, whose wave number is k — 0.002. A is the wavelength of this eigenmode, 
A = 2ir/k. The contour of the density and the temperature eigenfunctions are shown in (a) and 
(b), respectively, by grey scale, where the brighter (darker) region corresponds to the larger 
positive (negative) value. The arrows in (a) indicate the corresponding velocity eigenfunction. It 
shows that the grains flow from the brighter region into the darker density region, namely, the 
density perturbation decays. Note that the temperature perturbation is negative at the region 
of the positive density perturbation. 
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Figure 11. The eigenfunctions of the unstable mode for <j> = 0.05, 9 = 16°, and vo(0) = 0.05 
(PI in figure^, whose wave number is k = 0.002. The contour of the density eigenfunction and 
the velocity eigenfunction in (a) show that the grains flow from the region of negative density 
perturbation into that of positive perturbation, namely, the perturbation is amplified and results 
in the formation of the density wave. The temperature perturbation in (b) is negative at the 
region of positive density perturbation. 



in k for small k. As vq(0) is increased, the maximum value of Re [a] decreases and becomes 
negative for all k (figure |7{b) and (c)). 

The dispersion relations for the unstable modes at P4 (in the regime A) and at P5 (in 
the regime C) are shown in figure El for (9,v (0)) = (14°, 0.05) (a) and (17°, 0.10) (b), 
respectively. The instability occurs for the long-wavelength perturbation, and both of the 
dispersion relations show that a(0) = and the growth rate quadratic in k for small k, 
i.e., Re[a(fc)] oc k 2 ; these features are the same as those in the unstable regime B. 

The least stable eigenmodes for 6 = 16° at k = 0.002 are shown for the two cases: 
the stable case of ^o(0) = 0.10 (P3 in figure El in figure Hul and the unstable case 
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Figure 12. The stability diagram for h vs v with <f> = 0.10, where the unstable (stable) regimes 
are shown by grey (white) regions. We find three unstable regimes A, B, and C. The dashed line 
shows the boundary between the region of the non-monotonic density profiles and that of the 
monotonic density profiles (see also figure 2J. The dispersion relations at the points Q1~Q3 are 
shown in figure 

of ^o(0) = 0.05 in the regime B (PI in figure [SJ in figure ITT1 over one wavelength 
A = 2n/k. The contours in figures ITuT a) and lllf a) show the density eigenfunctions; the 
brighter (darker) regions indicate the positive (negative) regions, and the arrows represent 
the corresponding velocity eigenfunctions. In the both figures near y = 0, we see that 
the velocity perturbations point to the positive (negative) x direction in the regions of 
negative (positive) density perturbation. The contours for the corresponding temperature 
eigenfunctions plotted in figures HuT b) and lTTT b) show that the regions where the density 
perturbation is negative (darker regions in figures ITUT a) and II If a)) roughly correspond 
to the positive temperature perturbations (brighter regions in figures HUT b) and lllf b)). 

The difference between the stable mode (figure fTTTf) and the unstable mode (figure ITT)l 
is seen if we focus on the divergence of velocity perturbation. In the case of the stable 
mode for zaj(O) = 0.10 (figure ITuT a)). the grains flow into the region where the density 
perturbation is negative (see the region around x « A/2 and y ps 2), thus the density 
perturbation has negative feedback. On the other hand, in the case of unstable mode for 
vq (0) = 0.05 (figure HTT a)). the grains flow into the region where the density perturbation 
is positive (see the region around x « and y w 3); As a result, the perturbation grows 
and eventually causes the nonlinear density wave. 

The eigenfunctions of the unstable modes at P4 in the regime A and at P5 in the regime 
C show some different characteristics from those at PI in the regime B. Reflecting the 
difference in the steady flow between the regimes A and C and the regime B, the denser 
parts of the density eigenfunctions roughly correspond to the part where the velocity 
fluctuation has the positive component in the direction parallel to the mean flow. For all 
three regimes, however, the unstable modes show that the grains flow into the region of 
the positive density perturbation; this suggests that the instability leads to the density 
wave in the regimes A and C as in the regime B. 

4.2.2. The case of (f> = 0.10 

The stability diagram for <j> = 0.10 is shown in figure ITSl As in the case of <j> = 0.05, 
there are three unstable regimes, but the qualitative difference is that the unstable regime 
C at the large h and large v region contains a part of the region of the monotonic 
density profiles as well as that of the non-monotonic density profiles. The dispersion 



14 N. Mitarai and H. Nakanishi 

(a) (b) (c) 

Re(a) Re(a) Re(a) 




Figure 13. The dispersion relations of the least stable mode with <f> = 0.10 and 9 = 20° for 
iA)(0) = 0.10 (a) 0.15 (b), and 0.17 (c), which correspond to the points Ql, Q2, and Q3 in figure 
1121 respectively. 



relations around the boundary of the regime C are shown in figure ED f° r 6 = 20° and 
z^o(0) = 0.10(a), 0.15(b) and 0.17(c), which correspond to the points Ql, Q2, and Q3, 
respectively, in the stability diagram El It is seen that the instability occurs against the 
long-wavelength perturbation. 

The unstable eigenmodes for 9 = 20° and i/ (0) = 0.17 (Q3 in figure[I3) with k = 0.002 
are shown in figure El The density and the velocity eigenfunctions indicate that the 
grains flow into the the region of positive density perturbation (see the region around 
i » A/2, j/ w 0; the magnification is shown (c)), thus the density perturbation will grow. 
The difference from the case in figure El is that the velocity perturbation at the floor 
(y = 0) is in the positive x direction in the region of the positive density perturbation; 
namely, the particles flow faster in the region where they get dense. 



5. Discussions 

We have calculated the steady flow solutions and examined their linear stability under 
the longitudinal perturbation in the cases of <f> = 0.05 and = 0.10. The linear stability 
analysis revealed that there are three unstable regimes A, B, and C in the (h, v) plane 
in both of the cases. The regimes A and B are at the small v region, while the regime 
C is at the large v region. The difference between the regimes A and B is the density 
dependence of the flow velocity; the denser flow is faster in the regime A, while the flow 
with lower density is faster in the regime B. The regime C lies within the region where 
the density profile is non-monotonic in the case of <j) = 0.05, while the regime C includes 
a part of the region of the monotonic density profiles for <f> = 0.10, although the region 
of the regime C roughly corresponds with that of the non-monotonic density profiles. 
In all the regimes, the dispersion relation of the unstable mode shows the features that 
(i)ev(O) = and (ii)Re[a(fe)] oc k 2 for small k. The obtained unstable eigenmodes suggest 
that the instability causes the density wave. 

In this section, we give some discussions on the mechanism of the instability, com- 
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Figure 14. The eigenfunctions of the unstable mode for = 0.10, 6 — 20°, and ^o(O) = 0.17 
(Q3 in figure ^J, whose wave number is k = 0.002. The contour of the density eigenfunction 
and the velocity eigenfunction in (a) show that the grains flow into the the region of the positive 
density perturbation. In (b), the contour of the temperature perturbation is shown, (c) is the 
magnification of (a) around x ~ A/2, y « 1. 



parison with simulations, and the relationship with the transverse instability studied by 
Forterre k Pouliquen(2002]| 



5.1. Mechanism of the instability 

5.1.1. Kinematic wave and the long-wave instability 

The long-wavelength instability which results in the density wave is well known for 
quasi-one dimensional flows, such as the wave formation in the film flow (Smith 1993;Ooshida 
1999), the density wave formation in the granular flow in a narrow vertical pipe (Raafat, 
Hulin & Herrmann 1996; Moriyama et aZ.1998), and the jam formation in the traffic flow 
(Kerner & Konhauser 1993; Bando et aZ.1995; Mitarai & Nakanishi 2000a, 2000b). This 
instability is closely related to the continuity of "density" p (the thickness in the case of 
incompressible fluid, the density per unit length along the pipe for the granular flow, or 
the density of cars for the traffic flow): p obeys the equation of continuity in the form 

ggQM) | dQ(x,t) _ Q 
dt dx 

where Q is the flux. The flux Q may be expressed in terms of p and its spatial deriva- 
tives, but may also depend on time especially when the inertia effect exists. In the long- 
wavelength and the long-time limit, the effect of the spatial derivatives of p and the inertia 
on Q may be neglected. Then the flux is determined by the density; Q{x, t) = Qo(p(x, t)), 
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and equation (|5.1(l becomes 

ggjM)_ c ^.*) =0 with c= _dQo(p) (5.2) 
<9i dx dp v ' 

The wave that can be described by this equation is called the kinematic wave (Whitham 
1974). The effect of the inertia and spatial derivatives for the small but nonzero wave- 
length appears as the growth rate quadratic in k (Whitham 1974; Ooshida 1999), which 
causes the instability to yield the density wave. The actual form of the flux Q depends 
on a physical system, but we call the instability caused through this mechanism "the 



long-wave instability" after Smith(1993) Note that the unstable mode caused by the 
long- wave instability has the following three features for the complex growth rate a(k): 
(i) a(0) = 0, (ii) Re[a(fc)] oc k 2 for small k, and (iii) the phase velocity of the least stable 
mode c = Im[o!(fc)]/fc is given by — dQo/dp in the long-wavelength limit. 

The unstable modes obtained in the present analysis of the two-dimensional slope flow 
satisfy all of the features (i), (ii) and (iii); the features (i) and (ii) have been already 
pointed out in the text, and (iii) can be seen in figure OUb) in the case of <fi = 0.05, 
^o(O) = 0.05, and 9 = 16°. Actually, (iii) can be shown using (i) and the equation of 
continuity 1|2.1[) as shown in Appendix B. 

These features strongly suggest that the longitudinal instability in the slope flow is the 
long-wave instability of the kinematic wave in a quasi-one dimensional system for all of 
the regimes A, B, and C. 

One may think that, in the regime C, which lies in the region of the non-monotonic 
density profile in the case of (f> — 0.05 (figure EJ, the non-monotonic density profile 
might play a crucial role in the instability as in the case of the transverse instability of 
Rayleigh-Bernard type, where the convection occurs due to the non-monotonic density 
profile (Forterre & Pouliquen 2002; Carpen & Brady 2002); however, the fact that the 
regime C contains a part of the region of the monotonic density profile in the case of 
= 0.10 f figure IT^ suggests that the shape of the density profile does not determine the 
instability. It should be also noted that the present longitudinal instability occurs at the 
long wavelength while the transverse instability appears at the finite wavelength. 

5.1.2. One dimensional model 

In order to confirm that the longitudinal instability is the long-wave instability of 
the one-dimensional kinematic wave, we now try to reduce our two-dimensional model 
into the one-dimensional model that preserve the major features of the original model. 
In spite of crudeness of our procedure, we will see the obtained one-dimensional model 
shows roughly the same stability diagram for the long-wave instability of the kinematic 
wave. 

The way we obtain the one-dimensional model is to integrate the original equations 
in y-dircction from to oo; the idea is similar to Valance & Pennec(1998) where the 



one-dimensional model has been obtained for the flow in a vertical chute by integrating 
the equations across the chute width. 

We define the one-dimensional density p{x,t), the average velocity u(x,t), and the 
average temperature T(x, i) as 



p(x,t) = I u(x,y,t)dy, (5.3) 
Jo 

/>oo 

p(x,t)u(x,t) = / v(x,y,t)u(x,y,t)dy, (5.4) 
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p(x,t)f(x,t) = / u(x,y,t)T(x,y,t)dy. (5.5) 
Jo 

The one-dimensional equation of continuity obtained by integrating l|2.1|l is 

d t p + d x {pu) = 0, (5.6) 

where d t and d x represent d/dt and d/dx, respectively. <|5.6[) is in the form of l|5.1|l with 
the flux Q(x, t) given by Q = pu. 

By integrating the x-component of the equation of motion (|2.2|) , we obtain 



d t (pu)+d x / v(x,y,t)u(x,y,t) dy 
Jo 

= psmO + V yx (v(0),T(0),u(0))-dx [ I E xx (u,T,u)dy ) , (5.7) 



which determines the evolution of Q = pu. The second term in RHS is the shear stress 
at the floor and comes from the integration of d y Y, yx . 

To simplify (|5.7J) . we make the following approximations; (i) Replace the second term 
in LHS J Q v{x,y,t)u(x,y,t) 2 dy by pu 2 , (ii) Neglect the velocity in the y-direction, v, 
and use the dilute limit expressions for the the pressure and viscosities in T* xx , namely, 
we replace "E xx by vT — (4/3) j 2(0) \fTd x u, where the first term comes from the pressure 
and the second term comes from the dilute limit of the shear viscosity. The second 
viscosity £(z/, T) is neglected because it is a higher order quantity in v (see (|2.7(l 'l. (iii) 
Estimate the integration of the shear viscosity term by multiplying the integrand by 
the decay length of the density (T/cos#), and replacing T by T, u by u, i.e., replace 
(4/3)/ a (0) f~(VTd x u)dy by H(T) ■ (A/3)f 2 (0)(Vfd x u), where H(f) = T/cos9. 

Then we obtain the equation 



d t {pu) + d x (pu 2 ) =psm9 + i:y X (p(0),T(0),u(0))-d x (p^) + (A/3)f 2 (0)d x (H(T)VTd x u). 

(5.8) 

For the shear stress at the floor T, yx (v(0), T(0), u(0)), we have the boundary condition 
of momentum balance JOJ), i.e. B yx (v(0), T(0),<u(0)) = -r]*(v(0), T(0))u(0). 

To close equations (|5.6fl and l|5.8|l . we need the relation between (^(0), T(0), u(0)) 
and (p, T, u), and the equation for the average temperature T. We simply assume that 
u(Q) = u and T(0) = T, and we use the empirical relation between ^(0) and p for the 
steady flows; namely, we take z/(0) = F(p), with the form of F(p) determined from the 
steady solution obtained numerically for a fixed 6 and different values of ^o(O) by using 
z^o(0) = F(po) with po — J °° vo(y)dy. We further assume that f is also determined 
by the one-dimensional density p, rather than to use the integrated equation for the 
temperature. The form of T — T{p) is determined from the steady flows, namely, we 
assume T = (l/p ) J °° v (y)T (y)dy = f(p ). 

Now we finally obtain the one-dimensional model in the following form; 

d t p + d x (pu) = 0, (5.9) 
p [d t u + ud x u] = a(p) [U(p) -u]- d x {pf(p)) 

+(4/3)/ 2 (O)0 x (H(f (fi)jf{fidj\ , (5.10) 

where 

a(p) = V *(F(p),f(p)), U(p) = ^. (5.11) 

In this model, the steady solution is given by p — pa = const, and u = U(po), and 
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the flux of the steady solution is given by go = PoU(po). The equation l|5.10(l with H5.11JI 
shows that the velocity U(po) is determined by the balance between the acceleration by 
the gravity and the drag force from the floor. 

This model is almost the same form as the traffic flow model proposed by Kerner fc Konhauser( 1993)1 



but with the different form of the function U(p). In the traffic flow model, U{p) is often 
called the "optimal velocity function", which defines the density dependence of the car 
velocity and is usually a decreasing function. On the other hand, in the case of the granu- 
lar flow on a slope, U(p) depends on the inclination angle and the boundary condition at 
the floor, and can take various forms as can be seen in the plot of po vs Qo « qo = poll (po) 
f in figures Efb) andHJb). 

The linear stability analysis of the steady solution can be performed analytically (for 
the traffic flow model, see e.g. Wada & Hayakawa 1998). It is easy to show that the 
instability condition Re [a] > yields 



(a( Po )U'( Po )) 2 k 2 > +fi(p )kA (pof(p )yk 2 , (5.12) 

where p(po) = (4 / 3) f% (O)-ff(T(po ) ) ^/T(po) and the prime represents the differentiation 
by po. The k — > limit of H5.12fl gives the stability criterion, 

(pot/'M) 2 > (poT(po))'. (5.13) 

The explicit form of the long- wavelength expansion of the dispersion relation for the least 
stable mode is given by 

a(k) = -(p U(po))'tk + \(poU'(p )f - {poT(po))'] k 2 + 0(k 3 ), (5.14) 

a{p ) L J 

and we see that the phase velocity of this mode in the long-wavelength limit is given 
by c = —(poU(po))' = —dqo/dpo, which shows that this is the kinematic wave. The 
instability arises when the coefficient of the k 2 term becomes positive, which occurs if 
the change of the velocity against density fluctuation is too fast compared to the effect of 
the pressure that reduces the density fluctuation. Note that the criterion of the instability 
H5.13|l does not depend on the shear viscosity term p(po) because it only appears in the 
fourth order term in k (see H5.12l) h the approximation for that term does not crucially 
affect on the criterion. 

The stability diagrams for 4> = 0.05 and = 0.10 obtained from this one-dimensional 
model are shown in figures lT5l In spite of the crude approximations used in the derivation 
of the one-dimensional model, the stability diagrams are qualitatively similar to those of 
the original model. The similarity further indicates that the density wave formation can 
be understood by the long-wave instability in quasi-one dimensional systems, like in the 
film flow and the traffic flow. 

It should be noticed that, in this one-dimensional model, the effect of the parameters 
in the original model such as e p , 0, and <E> are included more or less implicitly in the 
functional forms of U{p) and T(p) (</> also appears explicitly in l|5.11(l through 77*), and 
the y-dependences of the variables also affect U(p) and T (p) through the integration. 
Any changes that affect U and T result in the changes of the unstable regions (see the 
criterion (|5.13() ). although the nature of the instability remains the same. 

Before concluding this subsection, let us make a few comments on the works by 



Wang, Jackson, & Sundaresan(1997) and Valance & Pennec(1998) on the stability anal 



f The flux qo in the one- dimensional model is not exactly the same as the flux Qo in the 
original two-dimensional flow, due to the approximations used to derive U(p). 




Figure 15. The stability diagram for (a)(f> — 0.05 and (b)cj> — 0.10 obtained from the 
one-dimensional model. The unstable (stable) regimes are shown by grey (white) regions. The 
dashed line shows the boundary between the region of the non-monotonic density profiles and 
that of the monotonic density profiles for the original model. 



ysis of granular flow in a vertical chute using hydrodynamic models of rapid granu- 
lar flow; Wang et a/.(1997) performed the linear stability analysis numerically as in the 
present work, and |Valance fc Pennec(1998) analyzed the density wave by deriving a one- 
dimcnsional model from hydrodynamic equations for rapid flow. 

In the analysis of Wang et a/.(1997)l a parameter region has been found where the 
steady flow is unstable against longitudinal long-wavelength perturbation to form the 
density wave (figure 9 and 10 in Wang et al. 1997). This instability might be also under- 
stood as the long- wave instability observed in the present analysis. On the other hand, 
they also found the instabilities for finite wavelength perturbations, which has not been 
observed here. 

The analysis by Valance fc Pennec(1998)| shows clear similarity of the instability in the 
chute flow to the one in the slope flow. The one-dimensional model they obtained (we 
call it the VP model) has the mathematical structure and physical mechanism similar 
to those of our one-dimensional model: the velocity of the steady solution is determined 
by the balance of the gravitational acceleration and the drag force from the wall, and 
shows the long-wave instability for which the criterion is determined by the change of 
the velocity against density fluctuation and the pressure term. 

These results suggest that the instabilities in the vertical chute flow and the slope flow 
are in the same class. 



5.2. Comparison with the molecular dynamics simulations 

We find some qualitative agreements between the present results and our previous sim- 
ulations (Mitarai & Nakanishi 2001) as follows. 

Our simulations were performed for a fixed inclination angle and a particular roughness 
of the slope with the periodic boundary condition imposed along the flow direction. 
Within the examined parameter region, the steady flow shows the monotonic density 
profile (Mitarai, Hayakawa, & Nakanishi 2002), and the flow with lower density has 
higher velocity. 

It has been demonstrated that the density wave appears only in the long system with 
low enough particle density. We have performed the simulations for several sets of the 
slope length L and the particle number N. In the case of the particle density Na/L w 1.0 
(single layer), the clear density wave is not formed for L = 250. 5cr and L — 501cr, whereas 



20 



N. Mitarai and H. Nakanishi 



the density wave appears for L = 1002(7. Upon changing the density with a fixed system 
length L — 501cr, the density wave is formed when Na/L « 0.75, while the steady flow 
is stable for the denser cases with Na/L w 1.0 and 2.0. 

These tendencies of the simulations agree with the behaviors around the unstable 
regime B of the present model on the three points: (i)The flow with lower density is 
faster, (ii)The flow with lower density is less stable, and (iii) The critical wavelength for 
instability is very long. 

Regarding (iii), the critical wave length A c = 2-zr /k c is much longer than the particle 
diameter; A c w 900<7 for 8 = 16° and ^o(0) = 0.05, for example. This seems to be 
comparable with our simulation results, where the critical slope length L c was between 
500c and 1000a for Na/L w 1.0. We do not understand yet how such a small wave 
number arises, but we suspect that it comes from the long mean free path in the large y 
region where the density is low, namely, the particles flying over the clusters for a long 
distance prevent the growth of clusters in a smaller length scale. 

Based on the observations (i), (ii), and (iii), the parameters that we have simulated 
happens to be in the regime B, but if the simulations are performed with different den- 
sities, inclination angles, and/or boundary conditions at the floor, the behavior that 
corresponds to the regime A or C may be found. 



5.3. Comparison with the stability analysis against the transverse perturbations 



Forterre & Pouliquen(2002) examined the linear stability of the granular flow on a wide 



slope against the perturbation transverse to the flowing direction, in order to understand 
the regular streak pattern along the flow direction observed in their experiments (Forterre 
& Pouliqen 2001). They mainly focused on the parameter region where the non-monotonic 
density profile (or what they call the "inverted density profile" ) is observed, because they 
expected that such a flow would be unstable to form the vortex rolls from the analogy 
to the Rayleigh-Bernard instability. They have shown that, the flow is unstable against 
the transverse perturbation in the large part of the parameter region where the inverted 
density profile is observed. The unstable mode shows the vortex-like pattern, and they 
concluded that the streaks observed in the experiments result from the rolls of vortices. 
They also found that the flow with the monotonic density profile becomes unstable for 
some parameters, but the detail has not been reported. 

One of the differences between this transverse Rayleigh-Bernard type instability and 
the longitudinal long-wave instability appears in the length scale of the instability. The 
longitudinal instability occurs against the long-wavelength perturbations in the long-time 
behavior as can be seen in the dispersion relations, while the transverse instability occurs 
at a finite length scale comparable with the vortex roll. Our analysis shows that there 
is a parameter region where both of the instability may occur, around the region of the 
non-monotonic density profiles (the regime C). It should be interesting to investigate 
how the two instabilities interfere or not by the full three dimensional analysis. 

In the large inclination angle, Forterre & Pouliquen(2002) also observed the square 
lattice pattern. This phenomenon cannot be understood by the simple superposition of 
the long-wave instability and the Rayleigh-Bernard type instability, because the length 
scale of the long-wave instability is much longer than the one of the lattice pattern. 



6. Summary 

The steady flows and their linear stability are analyzed for the granular flow on a slope 
using the hydrodynamic model with the constitutive relations derived from the kinetic 



Linear stability analysis of granular flow on a slope 



21 



theory of inelastic spheres. We mainly focused on the relatively low density region where 
the density decays monotonically. 

The stability diagram shows the three unstable regimes A, B, and C in the both cases 
of 4> — 0.05 and <fr = 0.10. Two of the unstable regimes A and B are in the lower density 
region, and the regime C is in the high density region. The difference between the regimes 
A and B is that, the denser flow is faster in the regime A at the small h side, while the 
flow with lower density is faster in the regime B at the large h side. The regime C is 
in the large h and large v region; it lies within the region of the non-monotonic density 
profile in the case of </> = 0.05, while it contains a part of the region of monotonic density 
profiles in the case of <\> — 0.10, although the region of the regime C roughly corresponds 
with that of the non-monotonic density profile. In all regimes, the instability occurs for 
the long-wavelength perturbations and results in the formation of density wave. It has 
been found that the behaviors around the unstable regime B agree with the tendencies 
in our previous simulations of density wave formations. 

The dispersion of the complex growth rate a(k) has the features that (i)a(0) = 0, 
(ii)Re[a(fc)] oc k 2 for small k, and (iii)Im[a(fc)]/fc — — dQo/dpo in the long- wavelength 
limit. These strongly suggest that the instability is the long- wave instability of the kine- 
matic waves, which is often found in quasi-one dimensional flows. This is different from 
the transverse instability studied by Forterre fc Pouliquen( 2002)1 where the flow is un- 



stable at the finite wave number. 

In order to confirm that the instability is the long-wave instability of the kinematic 
wave, we simplified the original equations rather heuristically into a one-dimensional 
model, and showed that the long-wave instability occurs in the derived one-dimensional 
model. The stability diagram obtained from the one-dimensional model corresponds qual- 
itatively to the one obtained from the original equations. 

N. M. is grateful to Ooshida Takeshi for informative discussions. This work was par- 
tially supported by Hosokawa powder technology foundation and Grant-in- Aid for JSPS 
fellows. 



Appendix A. Linearized equations 

In this appendix, the linearized governing equations for the longitudinal perturbations 
<|4.1I) - (I4.3|) and the boundary conditions are given (Forterre & Pouliquen 2002, see also 
Alam & Nott 1998). The superscript denotes that the quantities are for steady so- 
lution. The subscript v and T denotes the partial derivatives by the variables; pi = 
dp{v,T) / dv\ u=Uo .T=T - The subscript y indicates the total differential with respect to y, 
namely, p° = dp°\u(y),T(y))/dy = [dp{v,T)/dv\ v = VOtT =T Q ]vo, y + [dp{v,T) / dT\ u=U(hT=T() ]T Q 
vo,y = d^o(y)/dy, and so on. The expressions of p{v,T), /i(i/,T), £(^,T), k(v,T), and 
T(v,T) are given in l|2~7)l . and £(f,T) = £ - 2^/3. Then, by inserting l|P |l -l{0 |l into 
H2.1|I - (I2.3(I through (|3.1(l - (|3.3(l . we obtain the following expressions; 

[d t + u d x ]i>i + [vod^Ui + [u 0t y + ^o<3 y ]t>i = 0, (A 1) 

[sin0 - p° u d x + uo^yti® + uo^nl y + ^o^/^yM 
+ [-v d t ~ v u d x + (£° + 2fi°)d 2 x + [i y d y + ^°d 2 y } Ul 
+ [-vou , y + iiyd x + (£° + AicOdzdj,] 

+ [-p^d x + Uo, w Mt + Uo^ylJ-Ty + u a,ylJ-T9y]Tl = 0, (A 2) 

[- cos6» - ply - pldy + u Q ,yii v d x ]vi 
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+ [$d x + fd x d y + fj,°d x d y ]ui 

+ [-vodt - u u d x + e y 8 v + £°d 2 y + 2/jP y d y + 2p d 2 + moc^i 
+ [-Pry ~ Prdy + uo.^MtW - 0, (A 3) 

+ [-P°d x + 2p°u 0ty d y ]u 1 
3 

- 2 U o T o,y - P°d y + 2^ Q u ,yd x vi 

~v d t - l^uoflx + ^d 2 x + K° y d y + n°d 2 y 
+K TyTo, y + K° T T ,yy + K^T , y d y -T° T + T\ = 0. (A 4) 

The linearized boundary conditions at y = are obtained from 1)2. 9J1 and 1)2.10)1 as 

= 4u / 7 Vi + JVVL (A 5) 



dyTi = ~fL 



-d>ui - -$T n 



- ft 



(A 6) 



Appendix B. The long- wavelength limit of the phase velocity 

We have obtained the dispersion relation of the least stable mode, a = a(k), and 
numerically found that it satisfies 

a(0) = 0. (Bl) 
In this appendix, we show that the phase velocity of this mode, c = Im(a(fc))/fc, satisfies 

dQ 



dpo 



(B2) 



in the long- wavelength limit {k —> 0), which suggests the mode is the kinematic wave. 
Here, p and Q are the one-dimensional density and the flux of the steady solution, 
respectively, defined in 1)3.15)1 . 

Note that the derivative dQo/dpo is taken within a family of solutions for a fixed 
inclination angle as below. For a given inclination angle, the density at the floor vo(0) = (3 
is a continuous parameter to specify a steady solution. To express j3 dependence explicitly, 
we rewrite the steady solutions as 

v{x,y,z,t) = v (y;f}), (B3) 
u{x,y,z,t) = My; 0,0), (B4) 
T(x,y,z,t)=T (y;p). (B5) 

Then dQo/dpo is given by 

dQo _ dQ (/3)/d/3 _ J™[d(p (y; /3)u (y; p))/dp]dy 



dp dp (/3)/d/3 f™[dv (y;[3)/dj3}dy 



(B6) 

In order to show l)B 2)1 . let us first express c by eigenfunctions of the mode. By lineariza- 
tion of the equation of continuity 1)2.1)1 using 1)4.4)1 . we obtain the following expression: 



av{y) = -ik (v (y; (3)u(y) + u (y; 0)v(y)) - d y (y (y; (3)v{y)). (B 7) 
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Integrating (|B 7jl from y = to oo, we have 

/>oo />oo 

a t>(y)&y = -ikl [is {y; (3)u(y) + u (y; f3)£>(y)]dy. (B8) 
J o 

Here, v{y),u(y),v(y), and T(y) depend on the wavenumber /c, and we expands these 
functions with respect to k, i.e., 

Hv) = "o(y) +iki>i{y) + .... , = fi (y) + iku x {y) + (B9) 

=v {y)+ikvi (y) + ..., T(i/) =fo(y)+ifcf 1 ( W ) + ...., (BIO) 

where (vi(y),Ui(y),Vi(y),Ti(y)) do not contain fc. The long-wavelength expansion of the 
dispersion relation a — a(k) which satisfies a(0) = is obtained from (|B 8|l when 
Jo ^o(y)dt/ 7^ 0, and we have 

, w ., Ja° Wo(y;P)My) + Mv,0)My)]dy . n ,, 2 ^ 

a(/c) = -ife-y h 0(fc ). (Bll) 

Jo My)dy 

Namely, c is given by 

fo°° Mj/;/3)w (y) +My;/3)£o(y)]d2/ 
Jo ^olyjdy 

in the long- wavelength limit. 

Now all we have to do is to express Xq = (i> , uo, vq. To) with respect to the steady solu- 
tion X (y;[3) = {vo(y'i fl), u {y; (3), 0, T (y; (3)). Let us write H2.i p -Q2.3 fl and the boundary 
conditions in the matrix form: 

r)X 

B—=N(X). (B13) 

Here, X — (y,u,v,T), B is a constant matrix, and N is a nonlinear operator. A steady 
solution X (y;P) satisfies 

N(X (y;f3)) = 0, (B 14) 

therefore, from 14.41) . X = (v,u,v,T) satisfies 

dN(X)) 



dX 



\X=X o (y;0) 



, B ,X = a{k)BX, (B15) 



and expanding (|B 15(1 by k with a(0) = 0, we have 

dN(X)) 



dX 



\x=x (y;P) X " = ° ( B16 ) 



for the lowest order of k. On the other hand, differentiating l|B 14jl by (3, we obtain 

dN(X)) l dX (y;f3) 



ox lx = x o(y^ dp 



= 0. (B17) 



It is plausible that the mode of zero eigenvalue for k = does not degenerate, because 
the mass is the only one conserved quantity (the momentum is lost at the floor, and the 



energy is dissipated.). Thus, from i|B 16|l and i|B 17|l . we obtain Xo oc dXo(y;(3)/d[3, or 
more explicitly, 

{ dv o {y,0) duo(y;p) dT (y;/3) \ 
K uo, Uo , To) oc [—g^—, —dp—'°> -^T~ ) ■ ( B 18 ) 
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Using ljB12|l and (|B18|) . we get 

Jo" [My^)(du (y;(i)/dp) + u (y;f3)(dp (y;p)/d(3)}dy 
I~(dvo(3r,0)/dl3)dy 

and by comparing this with IjB 6|) . we have 

dQp/d/3 _ dQ 
dpo/d/3 dpo ' 

which is 



(B19) 



(B20) 



REFERENCES 

Ahn, H., Brennen, C. E. & Sabersky, R. H. 1992 Analysis of the fully developed chute flow 

of granular materials. J. Appl. Mech. 59, 109-119. 
Alam, M. & Nott, P. R. 1998 Stability of plane Couette flow of a granular material. J. Fluid 

Mech. 377, 99-136. 

Anderson, E., Bai, Z., Bischof, C, Blackford, S., Demmel, J., Dongarra, J., Du Croz, 
J., Greenbaum, A., Hammarling, S., McKenney, A. & Sorensen, D. 1999 LAPACK 
Users' Guide, 3rd edn. Philadelphia, PA: Society for Industrial and Applied Mathematics. 

Anderson, K. G. & Jackson, R. 1992 A comparison of the solutions of some proposed equa- 
tions of motion of granular materials for fully developed flow down inclined planes. J. Fluid 
Mech. 241, 145-168. 

Bando, M., Hasebe, K., Nakayama, A., Shibata, A., & Sugiyama, A. 1995 Dynamical 
model of traffic congestion and numerical simulation. Phys. Rev. E 51, 1035-1042. 

Boyd, J. P. 2001 Chebyshev and Foulier Spectral Methods (Second Edition). Dover. 

Campbell, C. S. 1990 Rapid granular flows. Ann. Rev. Fluid Mech. 22, 57-92. 

Canuto, C, Hussaini, M. Y., Quarteroni, A. & Zang, T. A. 1988 Spectral Methods in 
Fluid Dynamics. Springer Series in computational physics. 

Carpen, I. C. & Brady, J. F. 2002 Gravitational instability in suspension flow. J. Fluid Mech. 
472, 201-210. 

Cundall, P. A. & Strack, O. D. L. 1979 A discrete numerical model for granular assemblies. 

Geotechnique 29, 47-65. 
Forterre, Y. & Pouliquen, O. 2001 Longitudinal vortices in granular flows. Phys. Rev. Lett. 

86, 5886-5889. 

Forterre, Y. & Pouliquen, O. 2002 Stability analysis of rapid granular chute flows: formation 
of longitudinal vortices. J. Fluid Mech. 467, 361-387. 

Goldhirsh, I 2003 Rapid granular flows. Annu. Rev. Fluid Mech. 35 267-293. 

Gottilieb, D., Hussaini, M. Y. & Orsag, S. A. 1984 Theory and applications of spectral 
methods. In Spectral Methods for Partial Differential Equations (ed. R. G. Voigt, D. Got- 
tilieb & M. Y. Hussaini). SIAM. 

Jaeger, H. M., Nagel, S. R. & Behringer, R. P. 1996 Granular solids, liquids, and gases. 
Rev. Mod. Phys. 68, 1259-1273. 

Jenkins, J. T. 1992 Boundary conditions for rapid granular flow: flat, frictional walls. J. Appl. 
Mech. 59, 120-127. 

Jenkins, J. T. & Richman, M. W. 1986 Boundary conditions for plane flows of smooth, nearly 

elastic, circular disks. J. Fluid Mech. 171, 53-69. 
Jenkins, J. T. & Savage, S. B. 1983 A theory for the rapid flow of identical, smooth, nearly 

elastic, spherical particles. J. Fluid Mech. 130, 187-202. 
Johnson, P. C. & Jackson, R. 1987 Frictional-collisional constitutive relations for granular 

materials, with application to plane shearing. J. Fluid Mech. 176, 67-93. 
Johnson, P. C, Nott, P. & Jackson, R. 1990 Frictional-collisional equations of motion for 

particulate flows and their application to chutes . J. Fluid Mech. 210, 501-535. 
Kerner, B. S. & Konhauser, P. 1993 Cluster effect in initially homogeneous traffic flow. 

Phys. Rev. E 48, R2335-R2338. 



Linear stability analysis of granular flow on a slope 



25 



LOUGE, M. Y. & Keast, S. C. 2001 On dense granular flows down flat frictional inclines. Phys. 
Fluids 13, 1213-1233. 

Lun, C. K. K. & Savage, S. B. 1986 The effects of an impact velocity dependent coefficient 

of restitution on stresses developed by sheared granular materials. Acta Mech. 63, 15-44. 
Lun, C. K. K., Savage, S. B., Jeffrey, D. J. & Chepurniy, N. 1984 Kinetic theories for 

granular flow: inelastic particles in Couette flow and slightly inelastic particles in a general 

flowfield. J. Fluid Mech. 223, 223-256. 
Malik, M. R. 1990 Numerical methods for hypersonic boundary layer stability. J. Comput. 

Phys. 86, 376-413. 

Mayer, E. W. & Powell, K. G. 1992 Viscous and inviscid instabilities of a trailing vortex. 

J. Fluid Mech. 245, 91-114. 
Mitarai, N., Hayakawa, H. & Nakanishi, H. 2002 Collisional granular flow as a micropolar 

fluid. Phys. Rev. Lett. 88, 174301. 
Mitarai, N. & Nakanishi, H. 2000a Convective instability and structure formation in traffic 

flow. J. Phys. Soc. Jpn 69, 3752-3761. 
Mitarai, N. & Nakanishi, H. 20006 Spatiotemporal structure of traffic flow in a system with 

an open boundary. Phys. Rev. Lett. 85, 1766-1769. 
Mitarai, N. & Nakanishi, H. 2001 Instability of dilute granular flows on rough slope. J. Phys. 

Soc. Jpn. 70, 2809-2812. 
Mitarai, N. & Nakanishi, H. 2003 Hard-sphere limit of soft-sphere model for granular mate- 
rials: Stiffness dependence of steady granular flow. Phys. Rev. E 67, 021301. 
Moriyama, O., Kuroiwa, N., Matsushita, M. & Hayakawa, H. 1998 4/3 law of granular 

particles flowing through a vertical pipe. Phys. Rev. Lett. 80, 2833-2836. 
Ooshida, T. 1999 Surface equation of falling film flows with moderate Reynolds number and 

large but finite Weber number. Phys. Fluids 11, 3247-3269. 
Pouliquen, O. 1999 Scaling laws in granular flows down rough inclined planes. Phys. Fluids 

11, 542-545. 

Prasad, S. N., Pal, D. & Romkens, M. J. M. 2000 Wave formation on a shallow layer of 
flowing grains. J. Fluid Mech. 413, 89-110. 

Raafat, T., Hulin, J. P. & Herrmann, H. J. 1996 Density waves in dry granular media 
falling through a vertical pipe. Phys. Rev. E 53, 4345-4350. 

Richman, M. W. 1988 Boundary conditions based upon a modified Maxwellian velocity distri- 
bution for flows of identical, smooth, nearly elastic spheres. Acta Mechanica 75, 227-240. 

Savage, S. B. 1984 The mechanics of rapid granular flows. Adv. Appl. Mech. 24, 289-366. 

Smith, M.K. 199 The mechanism for the long-wave instability in thin liquid films. J. Fluid 
Mech. 217, 469-485. 

Valance, A. & Pennec, T. L. 1998 Nonlinear dynamics of density waves in granular flows 

through narrow vertical channels. Euro. Phys. J. B 5, 223-229. 
Wang, C, Jackson, R. & Sundaresan S. 1997 Instabilities of fully developed rapid flow of 

a granular material in a channel. J. Fluid Mech. 342, 179-197. 
Wada, S. & Hayakawa, H. 1998 Kink solution in a fluid model of traffic flow. J. Phys. Soc. 

Jpn. 67, 763-766. 

Whitham, G.B. 1974 Linear and nonlinear waves. John Wiley & Sons. 



